#ifndef PARAM_AIRFOIL_HICKSHENNE_HPP
#define PARAM_AIRFOIL_HICKSHENNE_HPP

#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
#include <cstring>
#include <iomanip>
#include <vector>
#include <utility>
#include "../../../include/VarArray.hpp"
constexpr double EPS = 1.0E-15;
const double PI_NUM = 4.0 * std::atan(1.0);

template <typename T>
struct point
{
    T x, y, z;
};

template <typename T>
struct vecMag
{ // 最后一个为数值大小
    T nx, ny, nz, mag;
};
// 条件为假时 打印信息并终止程序运行
inline void Assert(bool b, const char *mess)
{
    if (!b)
    {
        std::cerr << mess << std::endl;
        exit(EXIT_FAILURE);
    }
}
//==============================================================
template <typename T>
class Airfoil
{

  public:
    // Airfoil(){};
    /**
     * @brief Construct a new Airfoil object
     * 
     * @param grid 翼型网格文件 首行为标志位(1,up)(0,low) 接下来为坐标(必须有y)
     *              测试用
     */
    Airfoil(const std::string grid);
    /**
     * @brief Construct a new Airfoil object
     *   up_grid 上表面网格坐标
     * @param up_num 上表面坐标个数
     * @param low_num 下
     */
    Airfoil(const T (*up_grid)[3], int up_num, const T (*low_grid)[3], int low_num);

    // 下面行为 不允许拷贝, 不允许移动
    Airfoil(const Airfoil &) = delete;
    Airfoil &operator=(const Airfoil &) = delete;
    Airfoil(Airfoil &&) = delete;
    Airfoil &operator=(Airfoil &&) = delete;
    ~Airfoil()
    {
        if (own)
        {
            delete[] up_coord;
            delete[] low_coord;
        }
    }

    /**
     * @brief 输出坐标到流中 可以为文件或者屏幕
     * 
     * @param os 输出流
     */
    void print(std::ostream &os) const;
    /**
     * @brief 返回翼型坐标首指针
     *  你可以使用数组形式访问坐标
     * @return point<T>* 
     */
    point<T> *dataUp() const { return up_coord; }
    point<T> *dataLow() const { return low_coord; }
    // vecMag<T> *dataNormal() const { return normal; }
    int sizeUp() const { return up_size; }
    int sizeLow() const { return low_size; }
    point<T> getTail() const { return tail; }
    point<T> getHead() const { return head; }

  private:
    // 坐标首尾相接 首尾重合,算两个点
    // point<T> *coord{};
    point<T> *up_coord, *low_coord;
    // 翼梢和前顶点坐标 首尾坐标点不改变
    point<T> tail{}, head{};
    // point<T> *normal;
    // 坐标点数
    int up_size{}, low_size{};
    // 是否拥有资源, 从文件中读时分配了(动态)内存,
    // 使用数组初始化时只是引用 没有分配内存
    bool own;
    // 下面这个没用上
    // T mmax(T a, T b) const { return (b < a) ? a : b; }
    // 设置翼梢和顶点的坐标
    void setHeadTail();
};

template <typename T>
Airfoil<T>::Airfoil(const std::string grid)
{
    own = true;
    std::ifstream ins{grid};
    Assert(ins.good(), "Assert: cannot open airfoil grid file!");

    int up_tag{-1};
    ins >> up_tag;
    // 上表面
    if (up_tag == 1)
    {
        ins >> up_size;
        up_coord = new point<T>[up_size];
        for (int i = 0; i < up_size; ++i)
        {
            ins >> up_coord[i].x >> up_coord[i].y >> up_coord[i].z;
        }

        ins >> up_tag;
        ins >> low_size;
        low_coord = new point<T>[low_size];
        for (int i = 0; i < low_size; ++i)
        {
            ins >> low_coord[i].x >> low_coord[i].y >> low_coord[i].z;
        }
    }
    else
    {
        ins >> low_size;
        low_coord = new point<T>[low_size];
        for (int i = 0; i < low_size; ++i)
        {
            ins >> low_coord[i].x >> low_coord[i].y >> low_coord[i].z;
        }

        ins >> up_tag;
        ins >> up_size;
        up_coord = new point<T>[up_size];
        for (int i = 0; i < up_size; ++i)
        {
            ins >> up_coord[i].x >> up_coord[i].y >> up_coord[i].z;
        }
    }
    ins.close();
    setHeadTail();
}

template <typename T>
Airfoil<T>::Airfoil(const T (*up_grid)[3], int up_num, const T (*low_grid)[3], int low_num)
    : up_size{up_num}, low_size{low_num}, own{false}
{
    up_coord = (point<T> *)up_grid;
    low_coord = (point<T> *)low_grid;
    setHeadTail();
}

template <typename T>
void Airfoil<T>::setHeadTail()
{
    T tmp{-1e6}, mtmp{1e6};
    int ihead, itail;
    for (int i = 0; i < up_size; ++i)
    {
        // 找尾点
        if (tmp < up_coord[i].x)
        {
            tmp = up_coord[i].x;
            itail = i;
        }
        // 找首点
        if (up_coord[i].x < mtmp)
        {
            mtmp = up_coord[i].x;
            ihead = i;
        }
    }
    head = up_coord[ihead];
    tail = up_coord[itail];
}

/**
 * @brief 这是根据坐标之间的距离获得 首尾坐标
 * 
 * @tparam T 
 */
/**
template <typename T>
void Airfoil<T>::setHeadTail()
{
    // 获取尾部
    int index{0};
    T tmp{-1e6};
    for (int i = 0; i < sz; ++i)
    {
        if (tmp < coord[i].x)
        {
            tmp = coord[i].x;
            index = i;
        }
    }
    tail = coord[index];
    // 获得首部
    T distance{0.0};
    for (int i = 0; i < sz; ++i)
    {
        tmp = pow(coord[i].x - tail.x, 2) + pow(coord[i].y - tail.y, 2);
        if (distance < tmp)
        {
            distance = tmp;
            index = i;
        }
    }
    head = coord[index];
}
**/
template <typename T>
void Airfoil<T>::print(std::ostream &os) const
{
    os << std::scientific << std::setprecision(15);
    os << up_size << "\n";
    for (int i = 0; i < up_size; ++i)
    {
        os << up_coord[i].x << "\t"
           << up_coord[i].y << "\t"
           << up_coord[i].z << "\n";
    }
    os << low_size << "\n";
    for (int i = 0; i < low_size; ++i)
    {
        os << low_coord[i].x << "\t"
           << low_coord[i].y << "\t"
           << low_coord[i].z << "\n";
    }
}

//======================================================================
#if 1
class HicksHenne
{
  private:
    // 上下表面 公式中的 en   坐标的指数
    double *up_en, *low_en;
    // 不存, 只是引用  设计变量 delta
    const double *up_delt;
    const double *low_delt;
    // 公式 中 sin的指数  上下表面布置的凸点个数
    double Exp;
    int up_pos_num, low_pos_num;
    /**
     * @brief 检查是否越界
     * 
     * @param pos xn 从0到一 不是翼型的坐标 是bump 的位置
     * @param n 个数
     */
    void check(const double *pos, int n)
    {
        for (int i = 0; i < n; ++i)
        {
            Assert((pos[i] < 1.0) && (0.0 < pos[i]),
                   "Assert: (pos[i] < 1.0) && (0.0 < pos[i]) not matched in HicksHenne");
        }
    }
    // 计算 en
    void setEn(double *en, const double *pos, int n)
    {
        for (int i = 0; i < n; ++i)
        {
            en[i] = std::log10(0.5) / std::log10(pos[i]);
        }
    }

  public:
    /**
     * @brief Construct a new Hicks Henne object
     *   上下表面的 bump 可以不一样
     * @param udelt 上表面的 delta
     * @param up  上表面 bump 的位置
     * @param nu  个数
     * @param ldelt 下表面 delta 
     * @param low  下表面 bump 的位置
     * @param nl  
     * @param EXP  sin 的指数
     */
    HicksHenne(const double *udelt, const double *up, int nu,
               const double *ldelt, const double *low, int nl, double EXP = 3.0)
        : up_en{new double[nu]}, low_en{new double[nl]},
          up_delt{udelt}, low_delt{ldelt},
          up_pos_num{nu}, low_pos_num{nl}, Exp{EXP}
    {
        check(up, nu);
        check(low, nl);
        setEn(up_en, up, nu);
        setEn(low_en, low, nl);
        // memcpy(up_pos, up, sizeof(double) * nu);
        // memcpy(low_pos, low, sizeof(double) * nl);
    }
    /**
     * @brief Construct a new Hicks Henne object
     *  上下表面的 bump 位置和个数一样
     * @param udelt 上表面 delta 设计变量
     * @param ldelt 
     * @param up_low_pos  bump 的位置
     * @param num  bump 个数
     * @param EXP 
     */
    HicksHenne(const double *udelt, const double *ldelt,
               const double *up_low_pos, int num, double EXP = 3.0)
        : up_en{new double[num]}, low_en{new double[num]},
          up_delt{udelt}, low_delt{ldelt},
          up_pos_num{num}, low_pos_num{num}, Exp{EXP}
    {
        check(up_low_pos, num);
        setEn(up_en, up_low_pos, num);
        setEn(low_en, up_low_pos, num);
        // memcpy(up_pos, up_low_pos, sizeof(double) * num);
        // memcpy(low_pos, up_low_pos, sizeof(double) * num);
    }

    ~HicksHenne()
    {
        delete[] low_en;
        delete[] up_en;
    }

    /**
     * @brief 重载函数运算符
     * 
     * @param aerofoil  上面的类的对象
     * @return std::vector<double>  与所给翼型坐标对应的 y 改变量
     *     [0,up_size) 为上表面的 [up_size, up_size+low_size) 为下表面
     */
    std::vector<double> operator()(const Airfoil<double> &aerofoil, VarArray<double, 2> &deriMat)
    {
        const int up_size = aerofoil.sizeUp();
        const int low_size = aerofoil.sizeLow();
        std::vector<double> deltaY(up_size + low_size, 0.0);
        // std::cout << deltaY.size() << "test \n";

        const point<double> *up_coord = aerofoil.dataUp();
        const point<double> *low_coord = aerofoil.dataLow();

        const int coord_num = up_size + low_size;
        const int design_varNum = up_pos_num + low_pos_num;
        // 重复的首尾定点也算了
        deriMat.ReConstruct(coord_num, design_varNum);

        // const vecMag<double> *normal = aerofoil.dataNormal();
        const point<double> tail{aerofoil.getTail()};
        const point<double> head{aerofoil.getHead()};

        double x{0.0}, tmp{0.0};
        double deltaX{tail.x - head.x};
        // 上表面
        for (int i = 0; i < up_size; ++i)
        {
            x = up_coord[i].x;
            // 排除首尾点
            if (EPS < (x - head.x) && EPS < (tail.x - x))
            {
                // 归一化
                x = (x - head.x) / deltaX;
                for (int j = 0; j < up_pos_num; ++j)
                {
                    tmp = std::pow(std::sin(PI_NUM * std::pow(x, up_en[j])), Exp);
                    deltaY[i] += up_delt[j] * tmp;
                    deriMat(i, j) = tmp;
                }
                // 下表面对上表面没有影响
                for (int j = up_pos_num; j < design_varNum; ++j)
                    deriMat(i, j) = 0.0;
            }
            else
            {
                deltaY[i] = 0.0;
                for (int j = 0; j < design_varNum; ++j)
                    deriMat(i, j) = 0.0;
            }
            // std::cout << up_coord[i].x << "\t";
        }
        // 下表面
        for (int i = 0; i < low_size; ++i)
        {
            x = low_coord[i].x;
            // 排除首尾点
            if (EPS < (x - head.x) && EPS < (tail.x - x))
            {
                // 归一化
                x = (x - head.x) / deltaX;
                for (int j = 0; j < low_pos_num; ++j)
                {
                    tmp = std::pow(std::sin(PI_NUM * std::pow(x, low_en[j])), Exp);
                    deltaY[i + up_size] -= low_delt[j] * tmp;
                    deriMat(i + up_size, j + up_pos_num) = -tmp;
                }
                // 上表面对下表面没有影响
                for (int j = 0; j < up_pos_num; ++j)
                    deriMat(i + up_size, j) = 0.0;
            }
            else
            {
                // 首尾点
                deltaY[i + up_size] = 0.0;
                for (int j = 0; j < design_varNum; ++j)
                    deriMat(i + up_size, j) = 0.0;
            }
        }
        return std::move(deltaY);
    }
};

/**
 * @brief 见构造函数和子函数
 *    delta 要与 bump 位置一一对应
 * @tparam T 
 * @param up_num 
 * @param low_num 
 * @param udelt 
 * @param up 
 * @param nu 
 * @param ldelt 
 * @param low 
 * @param nl 
 * @param EXP 
 * @return std::vector<double> 
 *   deriMat 为返回的 y 对  delta 的导数
 */
template <typename T>
std::vector<double>
GetDeltaY_notSame(const T (*up_grid)[3], int up_num, const T (*low_grid)[3],
                  int low_num, const double *udelt,
                  const double *up, int nu,
                  const double *ldelt, const double *low,
                  int nl, VarArray<double, 2> &deriMat, double EXP = 3.0)
{
    Airfoil<T> aerofoil(up_grid, up_num, low_grid, low_num);
    HicksHenne param_foil(udelt, up, nu, ldelt, low, nl, EXP);
    return std::move(param_foil(aerofoil, deriMat));
}

template <typename T>
std::vector<double>
GetDeltaY_same(const T (*up_grid)[3], int up_num, const T (*low_grid)[3],
               int low_num, const double *udelt, const double *ldelt,
               const double *up_low_pos, int num,
               VarArray<double, 2> &deriMat, double EXP = 3.0)
{
    Airfoil<T> aerofoil(up_grid, up_num, low_grid, low_num);
    HicksHenne param_foil(udelt, ldelt, up_low_pos, num, EXP);
    return std::move(param_foil(aerofoil, deriMat));
}
#endif
//==============================================================================

// 计算导数
#endif